#ifndef RSIM_LU_H
#define RSIM_LU_H

/* -------------------------------------------------------------------------- *
 * File: LU.h                                                                 *
 * Authors: Kishor Bhalerao                                                   *
 * Email : kishor8dm@gmail.com                                                *
 * Contributors:                                                              *
 * Email:                                                                     *
 * copyright (c) 2010 Authors.                                                *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify it    *
 * under the terms of the GNU General Public License as published by the Free *
 * Software Foundation, either version 3 of the License, or any later version.*
 *                                                                            *
 * This program is distributed in the hope that it will be useful, but        *
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
 * for more details. (http://www.gnu.org/licenses/)                           *
 * -------------------------------------------------------------------------- */

#include "Matrix.h"
#include "Vector.h"
#include "LinAlgDecl.h"
#include <algorithm>
#include <cmath>

using namespace std;

namespace RSIM{


/** LU Decomposition.

Note: This file was obtained from Template Numerical Toolkit (TNT).
Details of TNT can be found on the website http://math.nist.gov/tnt/
This file is minimally modified to make it compatible with RSIM.
The originial file name (available for download on TNT website) 
is jama_lu.h.

<P>
For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
unit lower triangular matrix L, an n-by-n upper triangular matrix U,
and a permutation vector piv of length m so that A(piv,:) = L*U.
If m < n, then L is m-by-m and U is m-by-n.
<P>
The LU decompostion with pivoting always exists, even if the matrix is
singular, so the constructor will never fail.  The primary use of the
LU decomposition is in the solution of square systems of simultaneous
linear equations.  This will fail if isNonsingular() returns false.
*/

class LU{
        
        private:
                /** Array for internal storage of decomposition.  */
                Matrix  LU_;
                int m, n, pivsign; 
                Vector_<int> piv;
                
                Matrix permute_copy(const Matrix& A, const Vector_<int>& piv, int j0, int j1);
                Vector permute_copy(const Vector& A, const Vector_<int>& piv);	
        public:

                /** LU Decomposition
                @param  A   Rectangular matrix
                @return     LU Decomposition object to access L, U and piv.*/
                LU(const Matrix &A);

                /** Is the matrix nonsingular?
                @return     1 (true)  if upper triangular factor U (and hence A) 
                is nonsingular, 0 otherwise. */
                int isNonsingular();

                /** Return lower triangular factor
                @param     L_Mat is an empty matrix */
                void getL(Matrix& L_Mat);

                /** Return upper triangular factor
                @param     U_Mat portion of LU factorization. */
                void getU(Matrix& U_Mat); 

                /** Return pivot permutation vector
                @return     piv */
                const Vector_<int>& getPivot()const{return piv;}

                /** Compute determinant using LU factors.
                @return     determinant of A, or 0 if A is not square. */
                double det();

                /** Solve A*X = B
                @param  B   A Matrix with as many rows as A and any number of columns.
                @return     X so that L*U*X = B(piv,:), if B is nonconformant, returns
                0x0 (null) array. */
                Matrix solve(const Matrix &B);


                /** Solve A*x = b, where x and b are vectors of length equal	
                to the number of rows in A. 
                
                @param  b   a vector (Array1D> of length equal to the first dimension
                of A. 
                
                @return x a vector (Array1D> so that L*U*x = b(piv), if B is nonconformant,
                returns 0x0 (null) array. */
                Vector solve(const Vector &b);
}; // Class LU

} // namespace RSIM
#endif
